
if (!file.exists("output")) {
  dir.create("output")
}

# load data
source("HNC_metadata_tidy.R")
source("data_handling.R")
source("risk_factor_regressions.R")

sbsCOSMIC <- read.csv("Input_signature_attributions/Input_SBS_COSMIC_attributions.csv", stringsAsFactors = T, row.names = 1) %>% tibble::rownames_to_column("donor_id")
dbsCOSMIC <- read.csv("Input_signature_attributions/Input_DBS_COSMIC_attributions.csv", stringsAsFactors = T, row.names = 1) %>% tibble::rownames_to_column("donor_id")
idCOSMIC <- read.csv("Input_signature_attributions/Input_ID_COSMIC_attributions.csv", stringsAsFactors = T, row.names = 1) %>% tibble::rownames_to_column("donor_id")

additional_df <- list(sbsCOSMIC, dbsCOSMIC, idCOSMIC)

cosmic <- additional_df[[1]]
for (n in 2:length(additional_df)) {
  cosmic <- cosmic %>% left_join(additional_df[[n]], by = "donor_id")
}

# merge with cosmic sigantures and tidy epi variables
data <- data %>%
  left_join(cosmic, by = "donor_id") %>%
  mutate(
    country = factor(country, levels = c("Brazil", "Colombia", "Argentina", "Greece", "Italy", "Czech Republic", "Slovakia", "Romania")),
    argentina = ifelse(country == "Argentina", "yes", "no"),
    greece = ifelse(country == "Greece", "yes", "no"),
    italy = ifelse(country == "Italy", "yes", "no"),
    brazil = ifelse(country == "Brazil", "yes", "no"),
    czechRepublic = ifelse(country == "Czech Republic", "yes", "no"),
    romania = ifelse(country == "Romania", "yes", "no"),
    subsite = factor(subsite, levels = c("Larynx", "OC", "OPC", "Hypopharynx")),
    age_group = factor(age_group, levels = c("55-65", "0-45", "45-55", "65-75", "75+")),
    age = as.numeric(scale(age_diag))
  )

# Select signatures present in >n% of samples
# Dichotomize the signatures (Presence (signatures>0)="1",absence (signatures==0)="0")
# If signature is present in more than 75% cases, dichotomize by the median
data <- dichotomizeSigs(metadata = data, sigsn = cosmic, sigcutoff = 2)

# Dichotomized signatures to analyze
signatures <- grep("^(SBS|DBS|ID).*_cat$", names(data), value = T)

confounders <- c("age", "sex", "region", "subsite", "alcohol_ever", "tobacco_ever")

# Risk factors and clinical features to include in regression analysis
parameters <- c(
  "sex", "stage", "age",
  "OC", "OPC", "Larynx", "Hypopharynx",
  "alcohol_ever", "tobacco_ever", "tob_alc",
  "hpv_pos_opc",
  "bmi_cat", "oral_cat", "hotdrinks"
)

regression <- LRmodel(data, signatures, parameters, confounders, pval = 1, make_plot = F)

write.csv(regression[["regression_param"]], "output/Supplementary_Table_13.csv", row.names = F)

